Part 1 - Data importing and tidying

Starting up

As per last week, download and unzip today’s data into a new folder and open the file “tutorial_10.Rproj”. This will open up R and R studio and make sure that you don’t have to worry about filepaths.

Setting up R to run our analyses

Next, try running the following code by copying and pasting (Ctrl+C, Ctrl+V) it into your R session’s console.

# The following code functions ensure that the packages (like sub-programs that R uses to do things) are installed and attached
for(package in c('lme4', 'lmerTest', 'tidyverse', 'foreign', 'readr', 'sjstats')) {
    if(!require(package, character.only=T, quietly=T)) {
        install.packages(package)
        library(package, character.only=T)
    }
}

The data

The data for this week’s exercise comes is the same as last weeks - it is an organizational psychology dataset described in Klein et al (2000), and is stored as siop_merged.csv.

Variables

The file siop_merge.csv contains 750 employee-level observations nested within 50 work-groups. Ignoring the variable group ID (grpid), there are eight standardized variables, where a higher score indicates a higher level (e.g., higher pay or more negative leadership behaviours). These data were collected by individual survey of the employees, so measure individual level perceptions, except for the variable “physen”, which is the group level “physical work enviroment”.

Table 1. The variable names and variables included in siop_merge.csv.

The first thing we need to do is to read in the data. Copy and paste the following into your R console to load today’s data.

siop_merged <- read.csv("siop_merged.csv")



The data

The first variables we’ll be dealing with today are “Job satisfaction”, “Positive affect” and “physical work enviroment (group level)”. As always its a good idea to plot your data. This data has three variables, so you can either plot it with facets, or use a 3D plot. 3D plots are usually useless unless they are interactive, where they can be pretty and sometimes helpful for deveoping an understanding of you data. Click and drag this one to have a look at your data.

Figure 1. A 3D scatterplot of Job satisfaction, positive affect and group Physical enviroment, colour indicates group membership.



Part 2 - Random intercepts model with a predictor at the group level

Let’s fit a random intercept model with one individual level predictor - we’ll start with positive affect - and one group level predictor - physical work enviroment.

This model with random intercepts for group and a single level 1 predictor (positive affect) may be written:

\(JobSat_{ij} = \gamma_{00} + \gamma_{10} PosAff_{ij}+ \gamma_{01}PhysEn_{0j} + u_{0j} + \epsilon_{ij}\)

This means each person’s job satisfaction is equal to the overall intercept ( \(\gamma_{00}\) ), a regression term times their positive affect ( \(\gamma_{10}(PosAff_{ij})\) ), a regression term times the group level predictor ( \(\gamma_{01}PhysEn_{0j}\) ), a random effect for their group (\(u_{0j}\), like an error term for group), plus an error term (\(\epsilon_{ij}\)) (i.e., the residual, as the MLM will not perfectly predict their score).

PA_PE_model <- lmer("jobsat ~ 1 + posaff + physen + (1 | grpid)", data = siop_merged)
Which should be read as:
Terms Meaning
jobsat ~ job satisfaction is predicted by
1 + An intercept parameter which is the same for everyone
posaff a regression coefficient times each person’s positive affect
physen a regression coefficient times each group’s physical enviroment
(1 | grpid) A random intercept for each group

Remember summary(PA_PE_model) is useful in getting a summary of your output out!

summary(PA_PE_model)

This model says that the relationship between positive affect and job satisfaction is the same in every group, but that the baseline - the intercept - can vary by group, and that we can predict the group intercepts using the group’s score on physical enviroment. To visulise this, we can plot job satisfaction by positive affect and draw a regression line for each group. The model does not actually specify an intercept for each group, but rather assumes that the intercepts come from a normal distribution with a mean of the overall intercept and a variance as specified in the random effects output above. The output we might be interested is the variance of these groups’ intercepts as opposed to their individual values.

Figure 2

Figure 2. Job satisfaction by positive affect, lines show the estimated relationship by group, colour indicates group membership.


Part 3 - Random intercepts and random slopes

Now we will include a random slope for the association between positive affect and job satisfaction as well as a group level predictor for the intercept, but without other predictors for the slope. This means that positive affect becomes both a fixed and a random effect.

This model with random intercepts and a random slope for group, a single group level predictor (physical work enviroment), and a single level 1 predictor (positive affect) may be written:

\(JobSat_{ij} = \gamma_{00} + \gamma_{10} PosAff_{ij} + \gamma_{01} PhysEn_{0j} + u_{1j}PosAff_{ij} + u_{0j} + \epsilon_{ij}\)

This means each person’s job satisfaction is equal to the overall intercept ( \(\gamma_{00}\) ), a regression term times their positive affect \(\gamma_{10} (PosAff_{ij})\), a regression term times the group level predictor \(\gamma_{01}PhysEn_{0j}\), a random effect for positive affect ( \(u_{1j}PosAff_{ij}\) ), a random effect for their group (\(u_{0j}\), like an error term for group), plus an error term (\(\epsilon_{ij}\)) (i.e., the residual, as the MLM will not perfectly predict their score).

RS_model <- lmer("jobsat ~ 1 + posaff + physen + (1 + posaff | grpid)", data = siop_merged)
Which should be read as:
Terms Meaning
jobsat ~ job satisfaction is predicted by
1 + An intercept parameter which is the same for everyone
posaff a regression coefficient times each person’s positive affect
physen a regression coefficient times each group’s physical enviroment
(1 + posaff | grpid) A random intercept and a random slope for the relationship between positive affect and job satisfaction for each group

Remember summary(RS_model) is useful in getting a summary of your output out!

summary(RS_model)

That this model says that the relationship between positive affect and job satisfaction changes in different groups by some estimated ammount, and that we can predict the groups’ true mean scores using their group’s physical enviroment. To visulise this, we can plot job satisfaction by positive affect and draw a regression line for each group. The regression line slopes and intercepts change by group, because we now have random slopes and intercepts, allowing the relationship between positive affect and job satisfaction to change by group. Conceptually, this model is better thought of as assuming that the group slopes and intercepts come from a normal distribution some overall mean (i.e., the fixed effects for the intercepts and slopes) and an estimated variance. The output we’re often interested in above and beyond the fixed effects is the variance of these groups’ intercepts, and the variance of the groups slopes.

Figure 3

Figure 3. Job satisfaction by positive affect, lines show the estimated relationship by group, colour indicates group.


Part 4 - Random intercepts and random slopes with a group level predictor

Let’s fit a random intercept and slopes model with one individual level predictor - again positive affect - and one group level predictor - physical work enviroment, and one group level predictor of the random effects.

This model with random slopes and intercepts for group, a single level 1 predictor (positive affect), and an interaction between workload and physical enviroment may be written:

\(JobSat_{ij} = \gamma_{00} + \gamma_{10} PosAff_{ij}+ \gamma_{01}PhysEn_{0J} + u_{0j} + \gamma_{11}(PosAff_{ij}\times PhysEn_{0J}) + \epsilon_{ij}\)

This means each person’s job satisfaction is equal to the overall intercept ( \(\gamma_{00}\) ), a regression term times their positive affect \(\gamma_{10} (PosAff_{ij})\), a regression term times the group level predictor \(\gamma_{01}PhysEn + u_{0j}\), a random effect for their group (\(u_{0j}\), like an error term for group), plus an error term (\(\epsilon_{ij}\)) (i.e., the residual, as the MLM will not perfectly predict their score).

RS_model_2 <- lmer("jobsat ~ 1 + physen + posaff + posaff * physen + (1 + posaff | grpid)", data = siop_merged)
Which should be read as:
Terms Meaning
jobsat ~ job satisfaction is predicted by
1 + An intercept parameter which is the same for everyone
posaff a regression coefficient times each person’s positive affect
physen a regression coefficient times each group’s physical enviroment
physen:posaff a regression coefficient times the product of physical enviroment times positive affect
(1 + posaff | grpid) A random intercept and a random slope for the relationship between positive affect and job satisfaction for each group

Remember summary(PA_PE_model) is useful in getting a summary of your output out!

summary(RS_model_2)

That this model says that the relationship between positive affect and job satisfaction changes in different groups by some estimated ammount, and that we can predict the groups’ true mean scores using their group’s physical enviroment, and that the relationship between positive affect and job satisfaction is influenced by the physical enviroment.

To visulise this, we can again plot job satisfaction by positive affect and draw a regression line for each group. The regression line slopes and intercepts change by group, because we now have random slopes and intercepts, allowing the relationship between positive affect and job satisfaction to change by group.

Figure 4

Figure 4. Job satisfaction by positive affect, lines show the estimated relationship by group, colour indicates group.